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ABSTRACT 

Although quasi-periodic oscillations (QPOs) have been discovered in different X-ray sources, their origin is 
still a matter of debate. Analytical studies of hydrodynamic accretion disks have shown three types of trapped 
global modes with properties that appear to agree with the observations. However, these studies take only 
linear effects into account and do not address the issues of mode excitation and decay. Moreover, observations 
suggest that resonances between modes play a crucial role. A systematic, numerical study of this problem 
is therefore needed. In this paper, we use a pseudo-spectral algorithm to perform a parameter study of the 
inner regions of hydrodynamic disks. By assuming a-viscosity, we show that steady state solutions rarely 
exist. The inner edges of the disks oscillate and excite axisymmetric waves. In addition, the flows inside the 
inner edges are sometimes unstable to non-axisymmetric perturbations. One-armed, or even two-armed, spirals 
are developed, which provides a plausible explanation for the high-frequency QPOs observed from accreting 
black holes. When the Reynolds numbers are above certain critical values, the inner disks go through some 
transient turbulent states characterized by strong trailing spirals; while large-scale leading spirals developed 
in the outer disks. We compared our numerical results with standard thin disk oscillation models. Although 
the non-axisymmetric features have their analytical counterparts, more careful study is needed to explain the 
axisymmetric oscillations. 

Subject headings: accretion disks — hydrodynamics — instabilities — waves 



1. INTRODUCTION 

Quasi-periodic oscillations (QPOs) are strong coherent fea- 
tures in the power density spectra (PDS) that are found in dif- 
ferent low mass X-ray binaries (LMXBs , see the r eview arti- 
cle by van der Klis 2006 and Remillard & McClin tock 2006) . 
Their frequencies range from ~ 300 to ~ 1200 Hz, suggesting 
that they correspond to accretion flows very close to the cen- 
tral objects (van der Klis 2000, 2006). Understanding the ori- 
gin of QPOs will, therefore, lead to precise measurements of 
black hole spins as well as direct test s of Einstein's equi valent 
principle ( lPsaltisll2004 but also see lPsaltis et"ani2008h . Al- 
though several models have been proposed in the last decade 
to explain some aspects of the observations, the origin of 
QPOs is still a matter of debate. 

Local analytical studies of hydrodynamic accretion disks 
have shown thre e types of trapped global modes (see the 
review article by lKatol l2001l hereafter KOI, and references 
therein, or more recently Wagoner 2008). The coupling be- 
tween inertial oscillations and pressure variations lead to two 
different frequencies. The higher frequencies correspond to 
inertial acoustic p-modes and the lower ones correspond to 
gravit y ^-modes dPerez et al.l 119971 : lOrtega-Rodriguez et al.l 
I2OO2I) . Oscillatio ns in the vertical dire ction are called cor- 
rugation c-modes (Silber gleit et al.ll2001l) . Although these os- 
cillations have properties that appear to agree with the ob- 
servations, analytical studies can take only the linear effects 
into account and they do not address issues of mode exci- 
tation, propagation, and decay. Moreover, observations of 
QPOs in black holes strongly suggest that resonances be- 
tween mod es play a crucial role (A bramowicz & Kluzniak 
llpOl; Kluz niak & Abramowiczl200lb . This fact h as led to the 
develo pment of non-linear models (for example. lKatdl2004l 
l2008allbh . 

In principle, numerical simulations can be used to test these 
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models, including both linear and non-linear ones. One would 
like to carry out simulations for many dynamical time scales 
to obtain accurate PDS an d study their tempora l varia bil- 
ity. However, even with the IShakura & Svunyaevl (Il973h a- 
viscosity (which greatly simplifies the problem) the effec- 
tive/turbulent Reynolds numbers in thin accretion disks are 
still beyond the value that one could perform direct numerical 
simulations with standard numerical methods. 

We have developed earlier a two-dimensional, pseudo- 
spectral algorithm i n order to mo del thin accretion flows 
around black holes ( Chan et aU '2005). Spectral methods are 
high-order numerical methods that have a very small numer- 
ical dissipation. For two-dimensional flows, they are at least 
an order of magnitude more efficient then (low order) finite 
difference methods. Our algorithm solves the hydrodynamic 
equations in the disk plane (i.e., with r and 0) by using the thin 
disk approximation. In addition to the a-viscosity, high-order 
artificial viscosity is implemented by using a spectral filter, 
which affects only the large wavenumber modes and preserves 
correct shock properties (see Ma 1998a b, for details). We use 
a buffer zone to absorb outgoing waves at the outer boundary. 
For the inner boundary, because the flow is always super-sonic 
(toward the central object), no explicit boundary condition is 
needed. We apply this algorithm to study the very inner re- 
gions of two-dimensional, viscous, polytropic accretion disks 
around a black hole with zero spin. 

We find that steady state solutions do not exist for most 
combinations of the parameters. The inner edges of the disk 
often oscillate and excite axisymmetric waves with propa- 
gates all the way to the outer boundary. Other than that, non- 
axisymmetric spirals are sometimes developed in the simula- 
tions and provide a plausible explanation of the QPOs. De- 
pending on the values of the parameters, these spirals can be 
classified as inner one-armed spirals or inner two-armed spi- 
rals; in the cases with extreme Reynolds numbers, the results 
are steady global one-armed spirals. Although detailed study 
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of all these features is beyond the scope of this paper, we will 
make connections between the observed numerical phenom- 
ena and the standard analytical models. 

This paper is organized as following. The assumptions and 
the governing equations of the problem are presented in the 
next section. In |3] we derive an approximate steady state 
solution, which is used to initialize our simulations. We sum- 
marize the results from our simulations in ^ and describe 
our temporal analysis method in |5] Detailed discussions 
of different important features are presented in the subse- 
quent sections, namely, the axisymmetric rings in ^ non- 
axisymmetric spirals in ^ and the steady global spirals in 
^ We describe the limitations of our study in ^ Finally, we 
conclude with a discussion and suggest future research direc- 
tion in Ho] 

2. ASSUMPTIONS AND EQUATIONS 

We are interested in the dynamics of thin accretion disks. 
For simplicity, we neglect the motion in the vertical direc- 
tion. By assuming a polytropic equation of state, the hydro- 
dynamics is fully described by only three quantities, namely, 
the column density E and the two velocity components in the 
disk plane vv, v^. The governing equations are thus given by 
the continuity equation 

— + V-(Sv) = (1) 
and the (two-dimensional) Navier-Stokes equation 

S-^ + S(v-V)v = -VP+Vr+S]g. (2) 
ot 

In the above equations, P = KT,^ is the height-integrated ther- 
mal pressure with K and F being the polytropic constant and 
polytropic index. We denote by r the viscosity tensor, which 
(in Cartesian coordinates) is given by 



dvi ^ dvj \ 
dxj dxi J 3 



(3) 



We also define the sound speed = (FP/S)'/^ and the scale 
height H = Cs/fiK, where ilx is the "Keplerian" frequency 
defined below. 

We employ the lShakura & SvunvaevI a-prescription so that 
the effective kinematic viscosity is parametrized by the di- 
mensionless number a < 1 in the form t^ss = ac^H. The 
effects of general relativity is taken into account within 
the framework of the pseudo-Newtonian approximation, for 
which 

GM ^ 



(r-rs)'- 



(4) 



where rs = 2GM/c^ is the Schwarzschild radius and f is the 
unit vector in the positive r-direction. Therefore, the "Keple- 
rian" frequency is given by 



r!K: 



r-rs 



(5) 



We define our variables in terms of the natural length and time 
scales GM/c^ and GM/c^. The radial domain is chosen to be 
4 < r < 64 so it contains the innermost stable circular orbit 
(ISCO) at risco = 6. We run the simulations up to f = 20,000, 
which is about 325 orbital time scale at /"isco- All three dy- 
namic variables are output at every time unit, resulting into 
20,001 snapshots (including the initial condition at f = 0). 



There are three (physical) parameters in our equations, 
namely , the polytropic const ant K, the polytropic index F, 
and the lShakura & SvunvaevI a parameter. The accretion rate 
M, ins tea d of being an input parameter as in Milso m & Taarrj 
d 19961) or iMao et akl (l2008l) . is computed from our simula- 
tions. We choose three values for each of the parameters. 
Based on the thin disk assumption, we use K = 10~'^, 10"^^, 
and 10"^ with F = 1, 4/3, 5/3. We also take a = 0.1, lO'^ ^ 
1, which provide Reynolds numbers Re < 10^. The param- 
eter study, therefore, contains 27 different simulations (see 
Table[B. 

3. INITIAL CONDITIONS 

To avoid numerical instabilities at the inner boundary, we 
start the simulations at a state in which the flow close to the 
inner boundary is falling toward the central object. Because 
the numerical value of K we choose is small, we neglect the 
pressure and viscosity term, and solve an approximated so- 
lution of equations ((Tj) and Q. Using the superscript ^-^^ to 
denote the zeroth order quantities, the zeroth order hydrody- 
namic equations reduce to the very simple form 

^ ^-(rvf SW) = 0, (6) 



■ dr 



„(0) 



;i (0) 



GM 



' dr 



r (r-rs)^' 



dr r 
Integrating equation (|6]l, we obtain 

2rfE(''> = -M. 



(7) 



(8) 



(9) 



The accretion M here is nothing but a constant of integration. 

For the Navier-Stokes equation, there are two different 
classes of solutions. The first class is the Keplerian solution 
(KS): 



0, 



(10) 
(11) 



which is obtained by solving two algebraic equations; there 
is no constant of integration. The column density in this 
solution is completely arbitrary (as long as it is non-negative) 
and M = 0. 

The second class of solutions is the free-falling solution 
(FS): 



^,(0)2 

e= — 

2 

; (0) 



2r^' 



GM 

r-rs ' 



(12) 
(13) 



which are simply the conservation laws of energy and angular 
momentum. The constant of integration e has the meaning 
of specific energy and / has the meaning of specific angular 
momentum, respectively. The column density is obtained by 
combining equation (|9]l and ( fT2] i. 

For the region inside the ISCO, the rotational profile should 
follow FS. On the other hand, the disk is (almost) Keplerian 
far away from the central object. Hence, by requiring that 
the two classes of solutions match each other at some critical 
radius rc outside ISCO, the zeroth-order rotation profile is 



r < r^ 
r> rn 



(14) 
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It is convenient to define 
c/lnl^f") f 2 



(0) ^ . 



d\nr 



r < rr 



(3r-rs)/2{r-rs) , r>rc 



(15) 



as the background "shear". The value of is a free parameter 
in the zeroth-order solution, although it should not be too big 
compared to risco = 6GM/c^. 

Because of the thin disk assumption, the first order az- 
imuthal velocity is small v^'* compared to vj^*. Instead of us- 
ing the radial velocity equation, we can easily obtain the first 
order radial velocity vj.'' by integrating the angular moment 
equation, which gives 



(0): 



dr 



,(i)S(0)r!(")- 



(16) 



The constant of integration, L, describes the rate of angular 
momentum transport. Substituting for z^ss and recalhng the 
definitions of q and M, we obtain 

E(0)r . £ 
2TTqaTK-p^=M-^^. (17) 



Note that we differentiate il^, which originates from the disk 
scale height (and hence vertical gravity), from the zeroth order 
angular velocity il*^*'*. Taking the limit r oo and assuming 
that S decays slower than (the Shakura & Syunyaev solu- 
tion gives E ~ r"-'/'' for large r), the above equation reduces 
to 

(18) 



M = SnaTK lim 



We define the quantity C = lim^-^oo E*^*'^ ^^/r^K- In order to 
have a non-zero accretion rate, C should converge to some 
positive value. If the fluid is isothermal, this value simply sets 
the normalization of the density and does not affect the dy- 
namics. When r y 1, £ defines the sound speed and changes 
the behavior of the accretion disk. 

In standard accretion disk models, the value of L is usu- 
ally solved by a ssuming q = at some boundary layer (see 
iFrank et al] 120021) . However, for accreting black holes, it is 
more natural to assume that the density drop close to zero 
around the ISCO. Let e be a small parameter Equation ( fTTI l 
can be rewritten as 

.2 

L=(l-e)M- 



'ISCO 



GM 



nsco-rs V nsco 
The density outside the critical radius is then given by 
3 



1-(1 



C. 



(19) 



(20) 



nsco - 

The radial velocity can then be solved by the continuity equa- 
tion 

" (2.) 



„(1) : 



27rrE('') ' 

Inside the ISCO, gravity is the dominant effect. We solve 
the radial velocity by conservation of energy 

1 _ 1 

72 



ISCO 



-2GM 



1 



1 



-1/2 



(22) 



Here we denote by v,isco and hsco the radial velocity and 
specific angular momentum at the ISCO, respectively. The 



parameter e affects the above equation through Vrfsco- It 
should be chosen small enough so that the constraint vvisco <C 
hsco/nsco is always satisfied. 

The approximation described here introduces two extra pa- 
rameters. The first one is e, which controls the initial trans- 
port of angular moment. In our simulation, it is always taken 
e = 1%. The exact value of e is not important because the 
initial conditions evolve to the correct steady state solution in 
a few dynamical time scales. The second one is the limiting 
term C = lim,.^oo S*"* '"/JIk- It is chosen so E*"*(rinax) = 1 
and e = in equation ( l20l l. where r^ax = 64 is the outer radius 
of our computational domain. This makes E = 0(1) so the 
numerical values of K roughly represent c^. 

4. OVERVIEW OF RESULTS 

Figure [T] shows a representative set of snapshots from our 
parameter study. Each row contains four density contours of 
a particular simulation. The gray scale is linear, with white 
representing E = and black representing E = 6. From top 
to bottom, we can see the five most significant features: (/) A 
pure axisymmetric mode is excited near the ISCO and propa- 
gates outward. (//) An one-armed spiral develops only at the 
inner edge, while the other part of the computational domain 
is stable. (///) Axisymmetric rings appear and interact with 
the one-armed spiral at the inner edge, (iv) Two-armed spirals 
develop at the inner edge when the rest of the disk is stable. 
Finally, (v) the axisymmetric waves clean up the inner part 
of the disk. Inner spirals develop in this low density region 
and perturb the rings. The rings then break apart in late time 
and go through a turbulent transition state. After that, strong 
steady trailing spirals develop in the inner disks while global 
leading spirals develop in the outer disks. 

We will discuss these features in more detail in later sec- 
tions. Nevertheless, to provide a general picture, we review 
the standard models of thin disk oscillation. Following iKOlL 
we consider the small quantity h„,„(r) ^ P^'^^ /YP\ The sub- 
script m and n indicate different azimuthal and vertical modes, 
respectively. And f denotes the first order fluctuation in 
pressure. For simplicity, we will drop the subscripts m and n 
hereafter The small quantity h is governed by the equation 



'dr 



1 dh 

'■ — dr 



h = 0. 



(23) 



where Q = ui - mi^K an d k is the epicyclic frequency 
jNowak & Wagoneiill992l) . 

Assuming the perturbation h is local and has radial 
wavenumber k, the previous equation reduces to a dispersion 
relation 

{u^-K^){[d^-nn\) = [d^clk^, (24) 

where il± = c^/H = f^K is the vertical oscillation frequency. 
Setting « = for a two-dimensional disk, the p-mode is de- 
scribed by 

Co^ = {uj-mVLYf = K^ + clk^- (25) 



while the g-mode reduces to a trivial mode 



; w-mf^K = 0, 



(26) 



which describes pure rotations. We will argue in §6] that the 
axisymmetric rings seen in cases (;) and (;//) cannot be de- 
scribed by setting m = in equation ( l25T l. On the other hand, 
the non-axisymmetric spirals in iii), (Hi), and (iv) can be de- 
scribed by equation (|26] |. while the global spirals are the m=l 
p-modes. 
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Fig. 1 . — Density contours for a representative set of parameters (AT, F, a) at different times. Tlie gray scale is linear, with wliite representing S = and blaclc 
representing S = 6. From top to bottom, we have (i) pure axisymmetric waves, (('(') axisymmetric waves interacting with an one-armed spiral, (Hi) an one-armed 
spiral in the inner disk, (h>) two-armed spirals in the inner disk, and finally (i') a turbulent state in the inner disk, and one-armed spiral in the outer disk. See 
for an overview of these features. 
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TABLE 1 

Summary of results from the parameter study 



K 


r 


a 


Label 


Resolutions 






Init. carina 




Init. Wspii-ai 


'^.spiral 


Global Spiral 


10-3 


1 


1 


w 


129 X 32 


O.DUj 


8 390 


02 


015 












1-0.5 

i 


yp) 




9 790 
Z. / zu 


9 S9S 
Z. JZJ 


U.UJ 


098 
U.UZO 










1 

i 


U. i 




1 on V ^9 


U.oDU 


U. / D / 


decay 










10-3 


4/ J 


]^ 


(A\ 

K^) 


129 X 32 


1 1 470 


1 1 223 


01 


012 








in-3 


4/3 


1-0.5 
i 


(p-\ 
K^) 


iZV X jZ 


/^97 




~ U.UJ 


09 S 
u.uz J 








in-3 


4/1 
4/3 


n 1 

U. i 


if\ 
\J) 


1 9Q V "^9 
iZV X jZ 


1 1 Al 


1 019 
i .UJZ 


decay 










in-3 


S /I 
D/J 


1 
i 


((r\ 


1 9Q V ^9 

izy X jz 


1 /I Q^Q 


1 A OAQ 
i^+.Uty 


'-^ U.Ui 


OOQ 

u.uuy 








in-3 


S /I 


1-0.5 

i 


yM) 


1 9Q V ^9 
iZV X jZ 


^. J j'-t 


jUO 


~ n 09 
^ u.uz 


09S 
U.UZJ 








10-3 


S/l 
J/ J 


0.1 




129 X 32 


1.434 


1.297 


decay 










in-3-5 


1 

i 


1 

i 




9S7 V fiA 
Z J / X 04- 


9 790 

z. /zu 


9 (^"^S 
Z.OjJ 


--^ 09 
u.uz 


09 1 
U.UZl 










1 
1 


1-0.5 

1 




Z J / X 04- 


U.oDU 


U.o iU 


~ 01 

U.UJ 


01 1 
U.Uj 1 


--^0 10 oil 

~ U. lU— U. 1 1 


U.U70 






1 

1 


n 1 

U. 1 




9S7 V f\A 
Z J / X 04- 


A 979 
U.Z /Z 


n 9sn 

U.ZJU 


decay 






(1919 




10-3.5 


4/1 

4/ J 


]^ 




257 X 64 


3 627 


3 529 


P5 02 


018 








10-3.5 


4/1 
4/ J 


^-0.5 


in; 


257 X 64 


1 147 


i .U07 


Pi 03 


031 


~ f ^ OQ 19 
U.yJy — U. iZ 


?^ 01 




10-3.5 


4/1 
4/ J 


1 




257 X 64 


363 


338 






21 


206 




10-3.5 


5/3 






257 X 64 


4 534 


4 425 


P5 02 


018 








10-3.5 


5/3 


^-0.5 




257 X 64 


1 434 


1 367 


P5 03 


028 


~ 0(^ 15 


01 




10-3-5 


5/3 


0.1 


(r) 


257 X 64 


0.453 


0.425 


decay 




Si 0.20 


0.202 




10^ 


1 


1 


(s) 


513 X 128 


0.860 


0.842 


^ 0.02 


0.025 


decay 


weak 




10^ 


1 


1-0.5 


(t) 


513 X 128 


0.272 


0.264 


^ 0.03 


0.031 


decay 


weak 




10-4 


1 


0.1 


(u) 


513 X 128 


0.086 


0.082 


decay 






0.117 




lo-'t 


4/3 


1 


(V) 


513 X 128 


1.147 


1.049 


!^ 0.02 




decay 




steady 


10^ 


4/3 


1-0.5 


(w) 


513 X 128 


0.363 


0.349 


^ 0.03 




decay 




steady 


10^ 


4/3 


0.1 


(X) 


513 X 128 


0.115 


0.110 


decay 






0.114 




10-4 


5/3 


1 


(y) 


513 X 128 


1.434 


0.999 


!^ 0.02 




decay 




steady 


10^ 


5/3 


1-0.5 


(z) 


513 X 128 


0.453 


0.229 


^ 0.03 




decay 




steady 


10^ 


5/3 


0.1 


(@) 


513 X 128 


0.143 


0.138 


decay 




Si 0.10 


0.110 





In order to give a more quantitative description of the 
various simulations, we list some important numbers in Ta- 
ble[T] The first three columns are the parameters we explored. 
The fourth column gives a single letter label to each simu- 
lation. The fifth column lists the resolutions. Note that, as 
we decrease K, the kinematic viscosity j^ss decreases as well. 
Higher resolutions are then used to resolve the high Reynolds 
number flows ^ The sixth column lists the analytical accretion 
rate M^na calculated by equation ( fTsT i; while the seventh col- 
umn lists the numerical accretion rate Mnum computed from 
the simulations, 

Mnum =--fJJ ^"^r dt d(j}. (27) 

The time integrations are taken between t = 10000 and 20000 
so that T = 10000. The resulting Mnum are almost indepen- 
dent of the radius. The values listed in Table [T] are computed 
at r = 34, the mid-point of our radial domain, to minimize the 
boundary effects. The differences between Mana and Mnum are 
always less than 10% except for simulations (y) and (z). This 
agreement, together with the fact that the mean density pro- 
files are well described by our analytical approximations (see 
next section and Figure [Hi, imply that our initial conditions 
are reasonably close to the steady states. 

The eighth and ninth columns are the oscillation frequen- 
cies of the axisymmetric rings Wrfng- The values in the table 

' We have performed very high resolution benchmark simulations (not 
shown in this paper). These benchmarks are done in a smaller computa- 
tional domain with both constant kinematic viscosity and a-viscosity. Both 
the axisymmetric rings and inner spirals appear, implying they are robust fea- 
tures in two-dimensional viscous disks. However, for the cases with extreme 
Reynolds numbers, our computational resource does not allow us to ran a 
long enough super-high resolution simulation. Therefore, it is not clear to 
us whether the turbulent transient state is a numerical artifact or a physical 
phenomenon. See ^for detailed discussions. 



are computed at the ISCO although the frequencies are inpen- 
dent of the radius in all of our simulations. The eighth column 
lists rough descriptions of the axisymmetric rings in the early 
time of the simulations. If the axisymmetric rings decay away 
before t = 2,000, we identify them as relaxation effects of the 
initial conditions and mark "decay" in the table. If the oscil- 
lations remain strong at f = 2,000, we measure the time dif- 
ference between the density peaks and estimate the frequency 
by the equation cjiing = 27r/period. The values listed in the 
ninth column are more robust measurements of the frequen- 
cies. They are calculated based on (temporal) spectral anal- 
ysis of the simulations between t = 7,712 and 20,000. The 
details can be found in the next section. The symbol " — " in- 
dicates the power is distributed across a wide dynamic range, 
i.e., no specific frequency can be identified. Empty cell in- 
dicates the solution is stable so no significant frequency is 
found. 

The tenth and eleventh columns are the oscillation frequen- 
cies of non-axisymmetric inner spirals Wspii-ai- Note that, due 
to the interaction between inner spirals and the axisymmetric 
rings, the spirals do not have a unique frequency in simula- 
tions (k), (n), and (q). The differences between the peaks are 
uneven. We therefore write a frequency range for the esti- 
mated frequencies in the tenth column. The eleventh column 
lists the peak locations for the power spectra. Finally, the last 
column lists if a steady global spiral appears at late times in 
the simulations. 

5. TEMPORAL ANALYSIS 

In order to make the analysis more intuitive, we first plot 
the column density for each simulation at </> = as a contour 
of radius and time in Figure |2] The time axis is broken into 
three segments, namely, the early time [0, 1000], middle time 
[9500, 10500], and late time [19000,20000]. The gray scale 
is linear, with white representing S = and black representing 
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to,10"^ r=5/3,a = ,1 , (i 






r 1000 y'^OO ItWiO I "-O'Xt 2C'0Xt ') litCU ^f''Cr ti>5f-u I900i> I "HiOv r "i'JO S'l'tQ t v^oo J900"> DWj'"J 

Fig. 2. — Column density at </> = as contours of t (horizontal axis) and r (vertical axis). The gray scale is linear, with white representing and black 
representing 6. The horizontal axis is chopped into three segments, namely, the early time [0, 1000], middle time [9500, 10500], and late time [19000,20000]. 
The high density (dark) features that appear around r > 6 correspond to the axisymmetric rings. The weaker features that appear around r < 12, with fast temporal 
variability, are non-axisymmetric spirals. 
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Fig. 3. — The power spectra multiplied by lo, i.e., wP((aj), of the density as contours of frequency (horizontal axis) and radius (vertical axis). They are roughly 
the Fourier transforms of the panels in Figure [2] along the horizontal axis. The order of each panel is the same as in Figure[2] Because the first three azimuthal 
modes are the most important ones, we can use red, green, and blue to present m = 0, 1 , 2 modes, respectively. Deeper color represents higher power in the 
logarithmic scale. See 0and footnote 3 for the details of the coloring scheme. The dotted lines and dashed lines are the Keplerian and epicyclic frequencies, 
respectively. 
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S = 6. This is identical to the gray scale used in Figure[T] 

For all panels except (v), (w), (y), and (z), the high density 
(dark) features around r > 6 correspond to the axisymmetric 
rings. The slopes of these features, Sr/St, are their (radial) 
pattern speeds. The weaker features appearing around r < 
12 correspond to the non-axisymmetric spirals. Because the 
angular velocities are almost Keplerian in all simulations, it is 
easy to identify panels (k), (u), (x), and (@) with one-armed 
spirals, and panels (1), (o), and (r) with two-armed spirals, 
simply base on their temporal variability. 

Panels (v), (w), (y), and (z) are more complicated. The 
three time segments show that the properties of the flow are 
time-dependent. For simulations (v) and (w), the flows are 
laminar in the first two segments. The extra features that 
appear around r ~ 7 in the middle segment indicate that the 
axisymmetric rings interact with the very weak inner spirals. 
These extra features break apart and transit to a turbulent state. 
The same kind of transition happens in simulations (y) and 
(z). The almost time independent dark feature at r ^ 27 in 
panel (y) and r ^ 32 in panel (z) are the steady global one- 
armed spiral s (see the lowest row of Figure [T]!. 
We follow iKatol 's convention and write 



S(f,r,. 



E„,(w,r)e'('^'-"'"^> duj, 



where E„,(aj,r) is the Fourier transform/coefficient, 
power spectrum, defined by 

p„{uj,r) EE |i;„,(w,'")P, 



(28) 



The 



(29) 



measures the contribution of different frequencies for a par- 
ticular azimuthal mode as function of radius. Because the dy- 
namic variables are smooth^ and periodic along the azimuthal 
direction, the azimuthal modes are trivial to obtain by apply- 
ing discrete Fourier transform. However, the variables are not 
necessary periodic in time, or, at least the exact period is not 
known before doing the temporal analysis. The PDS obtained 
by squaring the discrete Fourier transform in a finite time do- 
main (i.e., periodogram) is not a good estimator of the actual 
power spectrum P^(uj,r). 

In order to impr ove the statistics, we use the standard "over- 
lap method" (see iPress et al.l [19921) to average over 11 peri- 
odograms. Each periodogram is computed over a section of 
2048 snapshots, with two subsequent sections overlapped by 
1024 frames. This requires 12 x 1024 = 12288 snapshots in 
total. We drop the first 7712 snapshots in our data to avoid 
analysing the initial relaxation. The Hann window function is 
applied to each segment before calculating the periodogram. 
The resulting variance in the PDS is approximately 11%. 

We plot the estimated ujP„,{LO,r) in Figure [3] The order of 
the panels are the same as Figure|2] Because only the first few 
azimuthal modes are important, we can use the red, green, and 
blue color channels to represent each of them. These three 
color channels are overlapped in a way that darker color rep- 
resents a higher power in a logarithmic scale^. For example, 
the red tone for panel (a), (d), and (g) indicates they are purely 

^ For most of the simulations, the Fourier transforms converge exponen- 
tially fast along the azimuthal direction. The exceptions are simulations (v), 
(w), (y), and (z), in which the flows go through turbulent transient states. Nev- 
ertheless, the spectral filters artificially cut off the high ni-modes and make 
the function smooth at the grid scale. 

^ Because each color channel is independent of each other, we are using 
the full three-dimensional color space. This is very different from the stan- 
dard one-dimensional color space (gray-scale, rainbow, etc). It is not obvious 
how to manipulate them by using standard vector graphics. Therefore, in Fig- 



m = axisymmetric modes. The brownish color for panel (j), 
(m), and (p) indicates they are dominated by m = modes, 
but have minor contributions from the m = 1 modes. The dark 
gray in panel (s), (v), and (y) indicates all m = 0, 1, and 2 
modes contribute. The sharp green and blue lines in panel 
(1), (o), (r), (u), (x), and (@) correspond to the fundamental 
modes and overtones of the one-armed and two-armed spirals, 
respectively. 

In addition, the dashed lines in Figure [3] are the epicycUc 
frequency K(r) and the dotted lines are the Keplerian fre- 
quency JIk in pseudo-Newtonian gravity. In contrast to the 
standard assumption, the oscillations in our simulations are 
not locally sub-Keplerian. Instead, the angular frequencies of 
both the axisymmetric rings and the non-axisy mmetric spi- 
rals are independent of the radius (see lMao et al. 2008). Their 
values are listed in the ninth and eleventh columns in Table [T] 

6. AXISYMMETRIC RINGS 

If the axisymmetric rings were only seen in our spec- 
tral algorithm, one would worry that they may be due to 
numerical artifacts of our code. The fact that these os- 
cillations were seen in previous studies of viscous spread- 
ing ri ngs (Papaloizou & Stanlev" '1986*; lOkuda & Mineshig^ 
119911: lOkuda et al. 1992; Godon 1995) as well as viscous 
disk models (Okudaetal. 1992; Chen & Taam 1992, 1995] 
IMilsom & Taam 1996; Mao et al. 2008; Reynolds & Mill3 
2008) suggest they are likely to be physical, at least under 
the thin disk assumption. Recentlv.l Revnolds & Milled (l2008h 
performed a detailed temporal analysis of both hydrodynamic 
and magnetohydrodynamic simulations. For their inviscid 
hydrodynamic simulations, similar axisymmetric oscillations 
are seen in the mid-plane, although the authors associate them 
with the g-modes. We are also aware of the current study of 
oscillations in tilted magnetohydrodynamic disks by Henisey 
& Blaes (private communication). 

Analytical studies of radial oscillations of axisym- 
metric modes i n accr etion disks have a long history. 
iBlumenthal et al.l (11984 ) studied the over s tability of axisym- 
metric oscillations and iLubow & Pringld (Il993h studied ax- 
isymmetric wave propagation in accretion disks. A detailed 
analytical study of the generation and propagation of these 
waves is beyond the scope of this paper. However, a sim- 
ple comparison sug gests that mo re physics is ne eded in the 
standard model, i.e., iKatol (1200 Ih . hereafter iKOlL in order to 
describe wave propagation in thin disks. 

Note that, when ojiing is well defined, it is always smaller 
than the maximum epicyclic frequency K^ax = '^('"peak) ~ 
0.0347, where rpeak = 2(2-1- \/3). This suggests the axisymmet- 
ric rings are generated as inertial waves and then propagate 
outward. Recalling the dispersion relation for axisymmetric 
inertial acoustic wave is up- = + c^k^, because > and 
the disk is thin, c^k^ ^ n\ ^ k^, it leads to a constraint on 
the local oscillation frequency u; ^ k. Therefore, KOI argue 



that a wave with frequency u), 



ring 



cannot propagate in 



the region [ri , r2], where ri and r2 are the Lindblad resonance 
points satisfying uj^„g = K{r). If the wave is excited far out 
and propagates inward, its wavelength becomes infinite at 
and is reflected back out. On the other hand, if the wave is ex- 

ure[3] we first "pixelize" the gray scale contour for each azimuthal mode. The 
red, green, and blue channel are then filled by the (maximum allowed) value 
255, resulting the overall red, green, and blue tones, for the m = 0, 1, and 2 
modes, respectively. To combine these different image, we use the graphic 
"multiply" operator. That is, each channel in the final image is a normalized 
pixel-wise product of the con'esponding channels of the three contours. 
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Fig. 4. — Typical density profiles at (/> = 0. The dotted line, dashed line, 
and dashed-dotted lines are profiles taken at t = 10000, 10100, and 10200, 
respectively. The thick solid line is the time average between t = 10000 and 
20000. The gray area marks the maximums and minimums of the waves. 

cited near the inner edge of the disk, it will be reflected back 
and forth between the inner edge and r\ . This phenomenon is 
called wave trapping. 

Taking simulation (j) as an example, coring = 0.021. The 
standard picture predicts no wave can propagate within the 
region [ri,r2] = [6.203, 12.587]; waves should be trapped be- 
tween the inner edge (near the ISCO) and r\. Panel (j) in 
Figure [3] clearly contradicts this prediction and violates the 
constraint w > k. Indeed, the power at the dominant mode 
increases in the forbidden region [ri,r2], and extends far to 
the outer region. KOI commented that the physical reasons 
for the difference between analytical models and numerical 
simulations are not clear. 

In Figure m we plot the density profile of simulation (j) at 
(/) = 0. The different dashed and dotted lines are profiles taken 
at t = 10000, 10100, and 10200. The thick soUd line is the time 
average between t = 10000 and 20000. The gray area marks 
the maximums and minimums of the waves. The amplitudes 
of the oscillations (for this particular simulation) are not small 
compared to the background. The peak is located around r = 
9.5, with peak densities about four times the mean density. 
Also, the wavelengths are comparable to the radius, i.e, k ^ 
\/r. 

To compare how the amplitudes vary with the parameters, 
we plot the quantity |max(I])/ (S) - 1 ^ ~ \h\^ in Figure|5] The 
solid, dotted, and dashed lines correspond to simulations (a), 
(j), and (s), respectively. The fluctuations in simulation (a) 
decays exponentially as a function of r. As the pressure de- 
crease, like in cases (j) and (s), the fluctuations become non- 
linear and the decay rates increase. 

There are two kinds of difficulties in applying the stan- 
dard thin disk oscillation modes. First, as shown in FigureHl 
the local linear approximation breaks down due to the non- 
linearity and large wavelength of the axisymmetric waves. 
Second, even in the cases where local linear approximation 
is valid [say simulation (a)]. Figure |5] shows the fluctuations 
are damped in large radii, hence the wavenumber k is com- 
plex. The standard approach, i.e., equation ( |23] |. does not take 
this into account. More physics perhaps related to viscosity 
should be included in the linear mode analysis in order to un- 
derstand the wave generation and propagation. 

7. NON-AXISYMMETRIC INNER SPIRALS 

Linear stability analysis of hydrod ynamic disks was an 
active area of study in the 1980's. iPapaloizou & Pringld 
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Fig. 5. — The function |max(S)/(S) — 1 1- ~ |/!pplotted in log scale taken 
from the period t = 10000 to 20000. The solid line is for simulation (a), the 
dotted line for (j), and the dashed line for (s). When the amplitude is small, 
i.e., for simulation (a), the fluctuation decays exponentially for r > 10. 

([1981 IT985h. i Goldreich & NaravanI (119851) and 
IPapaloizou & Pringle (1987) showed that inviscous tori 
are unstable to non-axisymme tric perturba ti ons. The global 
analysis was extended by |Blaesl (Il985l) . IGoldreich et al 
(1986), and iGoodman et all (Il987b : while iNaravan et al 
(119871) performed a local analysis based on the shear- 
ing sheet approxima t ion. Self- gravity was inclu ded in 
IGoodman & NaravanI ( Il988h . and IPapaloizou & LinI ( |1988|) 
showed that, if the effective viscous stress is a rapidly 
increasing function of the column density, non-axisymmetric 
perturbations grow because of viscous overstability. 

Table [T] shows that for one-armed spirals, o^spirai is al- 
ways very close to the Keplerian frequency at the ISCO, 
f^K('"isco) ~ 0.102; while for two-armed spirals cjspirai k, 
2f7K('"isco)- This observation strongly suggests that inner spi- 
rals are generated by comtating density fluctuations, which 
can be described by the trivial mode at the ISCO. The density 
fluctuations spiral inward as the SPBF suggests and create the 
rotating patterns in the inner disks. 

In Figure |6] we plot the contours of the very central re- 
gion of two simulations that have inner spirals but no signifi- 
cant axisymmetric rings. The gray scale is again linear, with 
white representing E = 0. Black represents S = 9 for the top 
row [simulation (k) with an one-armed spiral] and S = 6 for 
the bottom row [simulation (r) with two-armed spirals]. The 
time difference between each snapshot is 16, which is about 
a quarter of a period at the ISCO. Therefore, the contours at 
t = 19936 look almost identical to the contours at f = 20000. 

The solid lines in Figure |6] show the amplitudes of the 
column density E along the horizontal axis with the values 
matching the vertical axis. The amplitudes are parameter de- 
pendent. For example, the spirals have non-linear amplitudes 
for simulation (k); while for simulation (r), the amplitude is 
small, so perturbation methods are appli cable. This result 
is con sistent with the analytical study of IPapaloizou & LinI 
(Il988h . in the sense that the amplitude of the spirals is an 
increasing function of T. Also note that the solid lines are 
symmetric (because of the m = 2 mode) for the bottom row. 

8. TRANSIENT AND STEADY GLOBAL SPIRALS 

iKatol (Il983h showed that one-armed inertial-acoustic waves 
(with n = 0), which are just an eccentric deformatio n of the 
disk plane, have frequencies much lower tha n f^i^. lOgilvid 
unified a number of efforts (see reference in IOgilviell200ll) 
and derived a comprehensive set of evolutionary equations 
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Fig. 6. — Close up view of the central one-armed (top) and weak two-armed (bottom) spirals. The gray scale is linear in both rows with white representing 
S = 0. Black represents S = 9 for the top row and S = 6 for the bottom row. The time difference between each frame is 16, which is about a quarter of a period at 
the ISCO. Hence, the contours at f = 19936 is pretty much the same as f = 20000 because the flows rotate a full period. The solid Hnes show the column density 
along the horizontal axis, with the vertical axis matching the numerical values. The solid lines for the bottom row are symmetric because of the m = 2 modes. 



for eccentric disks. Later, S peith & Kiev! (l2003h performed 
a careful perturbation analysis'* and showed that the pertur- 
bations are independent of the dynamical time scales. These 
types of (almost) steady global spiral patterns are believed to 
be the cause of the V/R variations of Be stars (see iKOll and 

reference th erein). 

Following lKatol (Il983h . we consider the p-mode dispersion 
relation dZSl) with m = 1 , 



(30) 



For nearly Keplerian disks, k w JIk- If the disk is thin and the 
oscillation is global, i.e., Cjfc ^ Hk, the lower frequency can 
be approximated by 



1 



(31) 



Hence, for a simulation with a small enough sound speed, we 
expect our numerical solution to evolve to a state that agrees 
with the above analytical result. 

Form the bottom panels of Figure [T] we can briefly see the 
different stage of such an evolution. When t < 5000, axisym- 
metric waves are excited and clean up the inner part of the 
disk. Spiral structures can form in this very low density re- 
gion. The inner spirals interact with the axisymmetric rings. 
At f 10000, if the Reynolds numbers are high enough [sim- 
ulation (v), (w), (y), and (z)], both the axisymmetric rings and 
the inner spirals are able to break apart. The inner disks then 
go through a turbulent transition state. Although the turbu- 
lent transient states are potentially very important, we cannot 
draw strong conclusions from them. This is because these 
high Reynolds number flows are barely resolved even with 
513 X 128 g rids. Other concerns include the inner boundary 
effects (see iMcKinnev & Gammidl2002h and some numeri- 
cal issues are discussed in the next section. Nevertheless, the 

■* For a small kinematic viscosity, because the perturbation is in the high- 
est order term, the problem is singular. A stretching time transformation 
is needed to take in to account the two different time-scales, as done is 
ISpeith & KlevH2003l) . 



Strong trailing spirals that develop in the inner disks and the 
global leading spirals develop in the outer disks later on are 
very robust features. 

To focus on the trailing spirals, we plot the contours of the 
central region of simulation (x) and (z) in Figure |2l The con- 
tours are plotted using the same time scale as Figure|6]for easy 
comparison. The gray scale is again linear, with white repre- 
senting S = and black representing S = 18 for both rows. 
The solid lines show the column density along the horizontal 
axis, which show the spirals are very non-linear. Note that the 
spirals patterns are almost independent of time, as suggested 
in iKatQ (,1983i) and summarized in equation (|3TI ). 

9. LIMITATIONS 

Accretion disks involve different physics at different scales 
are very complicated systems. Although we have presented 
results based on careful numerical studies, one needs to be 
aware of the limitations in the simulations. In this section, 
we provide a list of some important concerns we have, and 
suggest several possible improvements for future studies. 

There are several numerical limitations that we need to be 
aware of. First, for the high Reynolds number cases, although 
513 X 128 grid points are used, we can barely resolve the ac- 
cretion flows. Our artificial viscosity may then have important 
effects on the flows, even though they are high order Second, 
in some of our simulations, the axisymmetric rings (together 
with the inflow, of course) clean up the inner region of the ac- 
cretion disks. The density is so low that a density floor is of- 
ten used to prevent negative values. Because this only happen 
in the inner region, it is likely that the artificially introduced 
mass would fall towards the central object very quickly. How- 
ever, it is not clear to us how important this effect is on the 
properties of our numerical solutions. Third, because we used 
an explicit Eulerian algorithm, the time steps are extremely 
small due to the fast flow and small grid at the inner bound- 
ary. Our high resolution simulations use more than a million 
time steps. The culmination of both truncation and round off 
errors may start polluting our numerical solutions. 

The above concerns suggest that, although the turbulent 
transient states found in our high Reynolds number Simula- 
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Fig. 7. — Close up view of the central strong one-armed spirals found in both simulation (x) [top row] and (z) [bottom row]. The contours are plotted at the 
same time as in Figure|6]for easy comparison. Note that the spirals patterns are almost independent of time, as suggested in Kato ( 1983). The gray scale is again 
linear, with white representing S = and black representing E = 18 for both rows. The solid lines show the column density along the horizontal axis, which show 
the spirals are very non-linear. 

tions [i.e., (v), (w), (y), and (z)] are potentially important, one 
should not draw strong conclusions from them. Better numer- 
ical algorithms with higher resolutions are needed to study 
them. 

In addition to the numerical issues, there are prob- 
lems associated with non-general relativistic hydrodynam- 
ics. iMcKinnev & Gammie (2002) performed a careful study 
of two-dimensional (in r-z plane) viscous hydrodynamic 
disks. They pointed out that numerical solutions in pseudo- 
Newtonian gravity depend on the location of the inner bound- 
ary. In our simulations, because the chosen sound speeds are 
small, the flows at the inner boundary are super-sonic. The 
boundary effect is therefore not important. Nevertheless, a 
general relativity version of our code will improve the numer- 
ical solutions and lead to more confident results. However, 
there is no unique way to formulate relativistic hydrodynam- 
ics with viscos ity. Methods like using flux-limited diffusion 
dWilson & Mathewsll2003]). or including an explicitly defined 
"finite time propagator" dKoide et al.ll2007h . are proposed. 

Finally, we have completely neglected the effect of mag- 
netic fields and adopted the a-viscosity prescription. It is well 
known that MHD turbulence in thick disk simulations tend 
to destroy cohere nt structures and wa sh away periodic signa- 
tures. Moreover, iPessah et al.l (l2008h pointed out that mag- 
netorotational instability (MRI) driven turbulence does not 
behave like a-viscosity. It is not clear whether equation ([T]i 
and ^ can capture the correct dynamics of accretion flows, 
even within a thin disk approximation. A study of oscilla- 
tions from MHD simulations of tilded disks [He nisey & Blaes 
(pri vate communication)] and thin MHD disks (IShafee et alJ 
l2008l) will therefore be very interesting. 

10. DISCUSSIONS 

In this paper, we have performed a parameter study of two- 
dimensional viscous accretion disks. We have found three 
types of robust features: the axisymmetric rings, the non- 
axisymmetric inner spirals, and the steady global one-armed 
spirals. When applying a temporal analysis of the column 
density, a large amount of power is contained in a few modes 
with frequencies that are independent of the radius. 



Although the properties of the axisymmetric rings do not 
completely agree with theoretical predictions, the frequen- 
cies coring are always smaller than the maximum epicyclic fre- 
quency Kmax- This Strongly suggests that they are axisymmet- 
ric p-modes. By taking a closer look at the simulations, we 
find out that the local, linear assumption breaks down in the 
inner disk. Moreover, based on the results of the simulations, 
we argue that the standard linear mode analysis is not enough 
to understand these rings. The fact that these inertial acoustic 
modes can propagate to large radii with constant frequencies 
have important implications — a frequency found in an ob- 
served oscillation can correspond to a large range of radii in 
the accretion disks. Therefore, one cannot use it to e stimate 
the size of the emission region (see iMao et al.ll2008l for re- 
lated discussions). 

The two types of non-axisymmetric modes have very dif- 
ferent properties. For the inner spirals, their frequencies are 
always very close to a small multiple of the Keplerian fre- 
quency at the ISCO, mili^iri^co)- Hence, they simply are ro- 
tating patterns originating at the ISCO and corresponding to 
the trivial g-mode J) = 0. The formation of the m = 1 and m = 2 
modes is probably due to viscous overstability. 

Although there are some situations in which both the ax- 
isymmetric rings and the non-axisymmetric spirals co-exist, 
the connection between these two features is not entirely clear. 
The excitation of rings usually cleans up the inner disks and 
prevents the inner spirals to develop. However, there are cases 
in which that the inner spirals perturb the axisymmetric rings. 
Steady global spirals are then excited. These steady patterns 
have properties that agree with the low frequency global m = 1 
p-modes. They are predicted as general features when the 
disks weakly deviated from Keplerian and have slow sound 
speed (KOI). The same feature can also be understood as ec- 
centric deformation. 

Modeling if QPOs can roughly be classified into two major 
schools. One conjectures that the observed oscillations cor- 
respond to the eigenfrequencies in the accretion disks. The 
other proposes that an inspiral stream of matter is responsible 
for the modulation. Our simulations show, depending on the 
parameter, both features can be excited in viscous disks. Al- 
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though the two frequencies are associated with axisymmetric 
rings and inner spirals are in ratio 1:3 instead of the observed 
2:3, our results provide new insights to understand and model 
QPOs. 

The author would Uke to thank Ramesh Narayan, Alexan- 
der Tchekhovskoy, Ruth Murray-Clay, and Robert Penna for 
useful discussions; Dimitrios Psaltis, Feryal Ozel, Martin Pes- 



Abramowicz, M. A., & Kluzniak, W. 2001, A&A, 374, L19 
Blaes, O. M. 1985, MNRAS, 216, 553 

Blumenthal, G. R., Lin, D. N. C, & Yang, L. T. 1984, ApJ, 287, 774 
Chan, C.-K., Psaltis, D., & Ozel, F. 2005, ApJ, 628, 353 
Chen, X., & Taam, R. E. 1992, MNRAS, 255, 51 
— . 1995, ApJ, 441,354 

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: 
Third Edition (Accretion Power in Astrophysics, by Juhan Frank and 
Andrew King and Derek Raine, pp. 398. ISBN 0521620538. Cambridge, 
UK: Cambridge University Press, February 2002.) 

Godon, P 1995, MNRAS, 274, 61 

Goldreich, P, Goodman, J., & Narayan, R. 1986, MNRAS, 221, 339 

Goldreich, P, & Narayan, R. 1985, MNRAS, 213, 7P 

Goodman, J., & Narayan, R. 1988, MNRAS, 231, 97 

Goodman, J., Narayan, R., & Goldreich, P 1987, MNRAS, 225, 695 

Kato, S. 1983, PASJ, 35, 249 

— . 2001, PASJ, 53, 1 

— . 2004, PASJ, 56, 905 

— . 2008a, PASJ, 60, 1 1 1 

— . 2008b, ArXiv e-prints, 808 

Kluzniak, W., & Abramowicz, M. A. 2001, Acta Physica Polonica B, 32, 
3605 

Koide, T., Denicol, G. S., Kodama, T., & Mota, P. 2007, Brazilian Journal of 

Physics, 37, 1047 
Lubow, S. H., & Pringle, J. E. 1993, ApJ, 409, 360 
Ma, H. 1998a, SIAM J. Numer. Anal., 35, 869 
— . 1998b, SIAM J. Numer. Anal., 35, 893 
Mao, S. A., Psaltis, D., & Milsom, J. A. 2008, ArXiv e-prints, 805 
McKinney, J. C, & Gammie, C. F. 2002, ApJ, 573, 728 
Milsom, J. A., & Taam, R. E. 1996, MNRAS, 283, 919 
Narayan, R., Goldreich, P., & Goodman, J. 1987, MNRAS, 228, 1 
Nowak, M. A., & Wagoner, R. V. 1992, ApJ, 393, 697 
Ogilvie, G. I. 2001, MNRAS, 325, 231 
Okuda, T., & Mineshige, S. 1991, MNRAS, 249, 684 
Okuda, T., Ono, K., Tabata, M., & Mineshige, S. 1992, MNRAS, 254, 427 



sah, Sukanya Chakrabarti for giving detailed comments on the 
manuscript; Gordon Ogilvie for pointing out important refer- 
ences on hydrodynamic overstability; Ken Henisey and Omer 
Blaes for kindly sharing their on-going work on disk oscilla- 
tions. The simulations presented in this paper were performed 
on a Beowulf cluster in the Physics Department at University 
of Arizona. The author is currently supported by an ITC fel- 
lowship. 



Ortega-Rodriguez, M., Silbergleit, A. S., & Wagoner, R. V. 2002, ApJ, 567, 
1043 

Papaloizou, J. C. B., & Lin, D. N. C. 1988, ApJ, 331, 838 
Papaloizou, J. C. B., & Pringle, J. E. 1984, MNRAS, 208, 721 
— . 1985, MNRAS, 213, 799 
— . 1987, MNRAS, 225, 267 

Papaloizou, J. C. B., & Stanley, G. Q. G. 1986, MNRAS, 220, 593 

Perez, C. A., Silbergleit, A. S., Wagoner, R. V., & Lehr, D. E. 1997, ApJ, 476, 

589 

Pessah, M. E., Chan, C.-K., & Psaltis, D. 2008, MNRAS, 383, 683 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. R 1992, 

Numerical recipes in C. The art of scientific computing (Cambridge: 

University Press, Icl992, 2nd ed.) 
Psaltis, D. 2004, in American Institute of Physics Conference Series, Vol. 

714, X-ray Timing 2003: Rossi and Beyond, ed. P. Kaaret, F. K. Lamb, & 

J. H. Swank, 29-35 
Psaltis, D., Perrodin, D., Dienes, K. R., & Mocioiu, 1. 2008, Physical Review 

Letters, 100, 091101 
Remillard, R. A., & McClintock, J. E. 2006, ARA&A, 44, 49 
Reynolds, C. S., & Miller, M. C. 2008, ArXiv e-prints 
Shafee, R., McKinney, J. C, Narayan, R., Tchekhovskoy, A., Gammie, C. F., 

& McClintock, J. E. 2008, ArXiv e-prints 
Shakura, N. I., & Syunyaev, R. A. 1973, A&A, 24, 337 
Silbergleit, A. S., Wagoner, R. V., & Ortega-Rodriguez, M. 2001, ApJ, 548, 

335 

Speith, R., & Kley, W. 2003, A&A, 399, 395 
van der Klis, M. 2000, ARA&A, 38, 717 

— . 2006, Rapid X-ray Variability (Compact stellar X-ray sources), 39-112 
Wagoner, R. V. 2008, New Astronomy Review, 51, 828 

Wilson, J. R., & Mathews, G. J. 2003, Relativistic Numerical Hydrodynamics 
(Relativistic Numerical Hydrodynamics, by James R. Wilson and Grant 
J. Mathews, pp. 232. ISBN 0521631556. Cambridge, UK: Cambridge 
University Press, December 2003.) 



